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Abstract PSRCHIVE is an open-source, object-oriented, scientific data analysis software 
library and application suite for pulsar astronomy. It implements an extensive range of 
general-purpose algorithms for use in data calibration and integration, statistical analysis 
and modeling, and visualisation. These are utilised by a variety of applications specialised 
for tasks such as pulsar timing, polarimetry, radio frequency interference mitigation, and 
pulse variability studies. This paper presents a general overview of PSRCHIVE function- 
ality with some focus on the integrated interfaces developed for the core applications. 



1 INTRODUCTION 



Within the pulsar astronomy community, a number of individuals and research groups have d evel- 
oped free ly-available software for a wide variety of purposes. The SlGPROcQ dLorimeri 12001 1) and 
PRESTCpl dRansom et al.L [20021 l2003h software packages are widely used in the search for new pulsars, 
DSPSpQenables real-time ph ase-coherent disp ersion removal (Ivan Straten & Bailesl 201 lb and compu- 



tation of the cyclic spe ctrum dDemoresll I20TTI) . both TEMPcFlt Taylo r & Weisbergll 19891) and TEMPO20 



dEdwards et all 120061) are used in th e analysis of pulse arriv al time estimates, and PSRCAlQ provides 
access to the ATNF Pulsar Catalog dManchester et all 120051) . This paper describes some of the basic 
and advanced functionality of the PSRCHlV^Zl project, which provides access to a comprehensive range 
of too ls commonly requir ed for the analysis of pulse profileQ data and the various metadata that describe 
them dHotan et 



nly requir 
all 12004 



Each of the above software projects has been refined through stages of early adoption and beta 
testing followed by more regular usage and feedback from the community. In turn, as the software 
grows more reliable and reputable, it reduces barriers to newcomers, thereby promoting growth in the 
discipline. These tools are now an indispensable resource to researchers in the field, and nearly all 
observational analyses of radio pulsar data published in the past few decades have relied on one or more 
of these packages. 



http://sigproc.sourceforge.net 

http://www.cv.nrao.edu/ sransom/presto 

http://dspsr.sourceforge.net 

http://tempo . sourceforge . net 

http://www.atnf.csiro.au/research/pulsar/tempo2 

http://www.atnf.csiro.au/research/pulsar/psrcat 

http://psrchive.sourceforge.net 

A pulse profile is any phase-resolved statistical quantity (e.g. flux density) integrated over one or more pulse periods. 
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Many of the most fundamental algorithms implemented by PSRCHIVE originate in the tim- 
ing analysis software d eveloped by the collaborators of the Parkes Southern Pulsar Survey (PSPS; 
[Manchester et all 1 1 996b . Over nearly two decades, these have been generalized, refined and incorpo- 
rated into a modular and extensible framework that employs object-oriented design principles and is 
primarily implemented using the C++ language^. To increase theportability of the code, it is currently 
managed using an open-source distributed version control systerrri and compiled using a cross-platform 
build systerro 

PSRCHIVE was developed in parallel with the PSRFlTSp] file for mat, which is fully compliant with 
the Flexible Image Transport Systen£3 (FITS: lHanisch et all [2Q0lh endorsed by NASA and the IAU 
and compatible with the recommendations of the International Virtual Observatory Alliance^ The 
modular, object-oriented design of the PSRCHIVE software separates the data analysis routines from 
file I/O, enabling the software to be easily extended to handle other data formats. In addition to main- 
taining backward-compatibility with the file format used by the original PSPS timing analysis software, 
PSRCHIVE currently provides s upport to read data i n eight different formats, including the European 
Pulsar Network flexible format (lLorimer et all Il998l ). the PRESTO PREPFOLD output format, and three 
different file formats used by pulsar instruments at the Arecibo Observatory. PSRCHIVE automatically 
determines the format of input data files, and all of the usage examples presented in this paper can be 
applied to any of the supported formats without modification. 

The portability and extensibility of PSRCHIVE fosters the incorporation of new features and func- 

j 1 1 — L — 1 1 1 1 1 

tionality, including rigorous polarimetric calibration (van Straten, 2004; Ord et al., 2004); various meth- 
ods of arrival time esti mation (e.g.lTkvlqrl.ll992ilHotan et all 1 20051: Ivan Stratenll2006l) : Faraday rotation 
measure determinati on (lHan et all 20061; iNoutsos et alll2008l) : propagation of the four th-order moments 



of the electric field (|yan Stratenl 12009 ); and statistical analysis of profile variability (iDemoresti 120071: 
lOslowski et all 12011 ). Development also continues on some of the more elementary algorithms, such 
as estimation of the off -pulse baseline and identification of the on-pulse region, and computation of the 
signal-to-noise ratio. 

The functionality of PSRCHIVE is distributed across a suite of specialised programs that are run from 
the command line in a typical UNIX shell environment. A subset of these programs, known as the Core 
Applications, provide access to general-purpose routines that are typically required for the majority of 
data analyses; these include 

- psredit - queries or modifies the metadata that describe the data set; 

- psrstat- derives statistical quantities from the data set and evaluates mathematical expressions; 

- psrsh - command language interpreter used to transform and reduce data sets; and 

- ps rplot - produces customized, diagnostic and publication-quality plots. 

The Core Applications employ a standard set of command line options and incorporate a command 
language interpreter that can evaluate mathematical and logical expressions, compute various statistical 
quantities, and execute a number of algorithms implemented by PSRCHIVE. Using the Core Applications 
and data that are available for download from the CSIRO Data Access Portal, this paper demonstrates a 
typical scientific workflow used to analyse observational data and produce pulse arrival time estimates 
for high-precision timing. As part of this demonstration, PSRCHIVE is used to perform radio frequency 
interference (RFI) excision, polarimetric calibration, and statistical bias correction. Throughout the pa- 
per, reference is made to the more extensive online documentation available at the PSRCHIVE web 
site. Detailed usage information is also output by each program via the -h command-line option. In 
sections [2] through [5] each Core Application is introduced with a description of the motivation and de- 
sign of the program followed by a demonstration of its use through a practical exercise. Sections [6] 



9 http://www.cplusplus.com 

10 http://git-scm.com 

1 1 http://sourceware.org/autobook 

12 http://www.atnf.csiro.au/research/pulsar/psrfits 

13 http://fits.gsfc.nasa.gov 

14 http://www.ivoa.net 
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through[8]demonstrate the use of PSRCHIVE to perform RFI mitigation, polarimetric calibration, arrival 
time estimation and bias correction. The concluding remarks in Section|9]include a description of some 
PSRCHIVE functionality that is currently under development and some ideas for future work. 



1.1 Observational Data 



The PSRCHIVE software processes observational data stored as a three-dimensional array of pulse pro- 
files; the axes are time (sub-integration), frequency (channel), and polarization (e.g. the four Stokes pa- 
rameters). The physical properties of the data are described by various attributes (also called metadata). 
A single data file containing one or more sub-integrations is typically called an archive. Sub-integration 
lengths may be as short as one pulse period (e.g. for single-pulse studies) or as long as desired (e.g. for 
the standard, or template profile, used for high-precision timing). 

The examples in this paper make use of nine days of observations of PSR J0437— 4715 made at 
20 cm with the Parke s 64 m radio telescope on 19 to 27 July 2003. Discovered in the Parkes 70-cm 
survey dJohnston et all Il993h . PSR J0437— 4715 remains the closest and brightest mill i secon d pulsar 
known; it has a spin period of ~ 5.7 ms, a pulse widt h of about 130 /is dNavarro et al and an 

average flux of 140 mJy at 20 cm dKramer et al l fl998l) . With a sharply-rising main peak and large flux 
density, it is an excellent target for high-precision pulsar timing studies. However, owing to the transition 
between orthogonally polarized modes of emission near the peak of the mean pulse profile, arrival 
time estimates der ived from observations of PSR J0437— 4 715 are particularly sensitive to instrumental 
calibration errors (ISandhu et ali Il997[ Ivan Stratenl t2006b . Furthermore, pulse -to-pulse fluctuati ons in 
the emission from this pulsar on timescales ranging from ~ 10 t o ~ 300 us (Jenetet all fl998l) place 
a fundamental limit on the timing precision that can be achieved foslowski et all 1201 1 ). These issues 
are discussed in more detail in Sections [7] and [8] which demonstrate the PSRCHIVE tools available for 
mitigating the impact of polarization calibration errors and correcting the bias due to self-noise. 

As described in the Appendix, these data are available for download from the CSIRO Data Access 
Portal using the Pulsar Search tool. A significantly reduced and more readily accessible form of the data 
is also available for download from Swinburne University of Technology. Throughout this paper, it is 
assumed that the full path to the directory containing the observational data downloaded from Swinburne 
is recorded using the $PSRCHIVE_DATA shell environment variable. 



2 QUERY AND MODIFY METADATA WITH PSREDIT 

The pulse profile data stored in a pulsar archive file are accompanied by metadata, or attributes, that 
describe various physical characteristics of the observation such as the source name, right ascension 
and declination, centre radio frequency and bandwidth of the instrument, start time and duration of the 
integration, etc. These attributes may be queried and modified using the psredit program, which is 
more fully documented onlinePl 

The keywords used by psredit to address the attributes in an archive are also understood and 
used by other Core Applications. For example, using psredit keywords, psrstat can perform 
variable substitution and evaluate mathematical expressions that include attribute values; similarly, 
psrplot can annotate plots (e.g. axes labels and titles) with attribute values as well as any of the 
mathematical expressions and statistical quantities provided by psrstat. This modularity of design 
allows the interfaces to the Core Applications to be remembered once and used often. 

Exercise: In the $PSRCHIVE_DATA/ mem directory, the receiver name is not set in any of the data files 
(* . ar). This can be verified by querying the receiver name attribute. 



cd $PSRCHIVE_DATA/mem 
psredit -c rcvr:name *.ar 



5 http://psrchive.sourceforge.net/manuals/psredit 
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Running psredit <f ilename> with no arguments will print a listing of every attribute in the file. 
Attribute names that end in an asterisk (e.g. int* : wt*) represent vector quantities. Specifying the 
attribute name without the asterisk will print a comma-separated list of every element in the vector; e.g. 
to print the centre frequency of every channel in every sub-integration 

psredit -c int : f req *.ar 



(These data files contain 1 sub-integration and 128 frequency channels.) To query the value of a single 
element, or range of elements, a simple array syntax can be used; e.g. 

psredit -c ' int : f req [ 34 , 5 6-60 ] ' *.ar 



will print only 6 values for each file. The single quotation marks in the above command are necessary 
to protect the square brackets from interpretation by the shell. 

Set the receiver name to MULT_1 using the standard output option^} to overwrite the original files. 

cd $PSRCHIVE_DATA/mem 

psredit -c rcvr : name=MULT_l -m *.ar 



The data files in $PSRCHIVE_DATA/mem/ are a mixture of three different types. Create two sub- 



directories called pulsar/ and cal/ then use psredit to query the type 
sort the files into the two sub-directories. 


attribute and use this to 


cd $PSRCHIVE_DATA/mem 




mkdir pulsar/ 




mkdir cal/ 




mv "psredit -c type *.ar | grep Pulsar | awk '{print 


$1 } ' ' pulsar/ 


mv * . ar cal/ 





The last line of the above commands places all three calibrator file types in the cal / sub-directory. 



3 EVALUATE DATA WITH PSRSTAT 

In addition to accessing the physical attributes that describe the observation, it is also useful to compute 
derived quantities, such as statistical measures that describe the quality of the data, the effective width 
of the pulse, the degree of polarisation, etc. A wide variety of derived quantities can be computed using 
the psrstat program, which is more fully documented onlinePI. Running psrstat <f ilename> 
without any command-line arguments will print a listing of every available quantity. 

Any of the quantities (attributes or computed values) provided by the psrstat interface can be 
substituted into mathematical expressions that can be evaluated (expressions to be evaluated are enclosed 
in braces). For example, to search for significant peaks in a series of single-pulse archives, query the 
maximum amplitude in all phase bins normalized by the off -pulse standard deviation, 

psrstat -c ' { $all : max/ $of f : rms } ' <filenames> 



To query the effective pulse width in microseconds, 

psrstat -c '{ $weff* $int [ ]: period* le6 } ' <filename> 



http://psrchive.sourceforge.net/manuals/guide/design/options.shtml 
http://psrchive.sourceforge.net/manuals/psrstat 
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psrstat -c snr pulsar/*. ar 



Note that, by default, psrstat computes the S/N of the profile in the first sub-integration, frequency 
channel and polarization (indexed by subint, chan and pol, respectively). Take one file and print 
the S/N in each frequency channel using the command line option to loop over an index. The number 
of frequency channels in the file can be queried with the nchan attribute. 

psredit -c nchan pulsar/n2003200180804 . ar 

psrstat -1 chan=0-127 -c snr pulsar/n2003200180804 . ar 



The psrstat program can be used in combination with other common tools (such as GNUPLOlP^I) 
to investigate problems and/or verify the quality of data. When doing so, it is practical to use the -Q 
command line option to print only the value of each attribute, instead of key=value. 

4 TRANSFORM AND REDUCE DATA WITH THE PSRSH INTERPRETER 

PSRCHIVE includes a command language interpreter that provides access to a large number of common 
data processing algorithms, including radio frequency interference mitigation and polarimetric cali- 
bration. Access to this interpreter is provided by the psrsh program, which may be used either as an 
interactive shell environment or as a shell script command processor as more fully documented online^ 
The psrsh interpreter is also embedded in the command line interfaces of the Core Applications 
(psrstat, psrplot, and psradd). These applications use the interpreter to execute preprocessing 
tasks on input data files as described in Section 2.3 Job preprocessor of the online documentation. 
The full list of available commands is listed by running psrsh -H. The first column of the output is 
the command name, the second column is a single-letter short-cut key in square brackets, and the third 
column is a short description of each command. 

Exercise: The S/N values output by psrstat in the previous section are those of only one polariza- 
tion, not the total intensity. The four polarization parameters stored in these files describe the elements 
of the coherency matrix: AA, BB, K[AB], and 3 [AB]. The total intensity, I = AA + BB is formed by 
the pscrunch command. Use the preprocessing capability of psrstat to print the S/N of the total 
intensity as a function of frequency for one file. 

psrstat -j pscrunch -1 chan=0-127 -c snr pulsar/n2003200180804 . ar 



Add the f scrunch command to print the S/N of the total intensity integrated across the entire ob- 
serving bandwidth for each file. 

psrstat -j pscrunch, f scrunch -c snr pulsar/*. ar 



Using single-letter short-cut keys, the above line is equivalent to 

psrstat -j pF -c snr pulsar/*. ar 



The output of the above command (combined with the -Q command line option) can be redirected to 
a file and GNUPLOT can be used to plot the variation of S/N as a function of time due to interstellar 
scintillation. 



http://gnuplot.info 

http://psrchive.sourceforge.net/manuals/psrsh 
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5 DISPLAY DATA WITH PSRPLOT 

Using the PGPLOtFI graphics subroutine library, psrplot produces both diagnostic and publication- 
quality plots as documented onlind^l The psrplot program is a highly configurable plotting tool for 
use during all stages of data analysis and manuscript preparation. To see a full listing of available plot 
types, use the -P command-line option. Each plot can be configured using a wide range of options, 
including selection of the range of data to be plotted (zooming), specification of plotting attributes such 
as character size and line width, and definition of plot labels. Run psrplot -A <name> to list the 
generic options that are common to most plots, and psrplot -C <name> to list the options that are 
specific to the named plot. 

Plot labels may include any of the attributes accessible via psredit and/or quantities and mathe- 
matical expressions computed by psrstat. For example, to produce a publication-quality plot of the 
total intensity profile with the S/N printed inside the top-right corner of the plot frame, 

psrplot <filename> -pD -jFp -c below : r=' S/N : $snr' -c set=pub 



Note that filename(s) need not necessarily be the last argument(s) on the command line. 

Exercise: Use psrplot to plot the phase-vs-frequency image of the total intensity of the pulsar signal 
in each file. 

psrplot -p freq -j p pulsar/*. ar 



By default, the dispersion delays between frequency channels are not corrected. Using the command line 
option to loop over an index, plot the phase-vs-frequency images of 5R[A_B] and Ss[AB] (pol=2 , 3), 
which effectively correspond to Stokes U and V. 

psrplot -p freq -1 pol=2,3 pulsar/*. ar 



These quantities vary with frequency due to an instrumental effect; the two orthogonal polarizations 
propagate through different signal paths with slightly different lengths, introducing a phase delay that 
varies linearly with radio frequency. 

Use the loop-over-index option to plot the total intensity profile (the plot type named flux or its short- 
cut D) as a function of frequency in one pulsar data file. 

psrplot -pD -1 chan=0- -jp pulsar/n2003200180804 . ar 



Note that chan=0- specifies the entire range without having to know the index of the last frequency 
channel. The pulse profile is significantly distorted at the edges of the band due to quantization error 
(also called "scattered power") that arises during analog-to-digital conversion using 2 bits/sample. The 
most severely affected channels will be excised in Section [6] 

Plot the Stokes parameters integrated over all frequency channels for each file. 

psrplot -p stokes -jF pulsar/*. ar 



The white line is the total intensity; red, green, and blue correspond to Stokes Q, U, and V, respectively. 
The Stokes parameters vary with time owing to the rotation of the receiver feed with respect to the sky 
(the parallactic angle). In Section|7] this effect will be exploited to model the polarization cross-coupling 
in the instrumental response and calibrate the data. 



http://www.astro.caltech.edu/ tjp/pgplot 
http://psrchive.sourceforge.net/manuals/psrplot 
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6 RADIO FREQUENCY INTERFERENCE AND INVALID DATA EXCISION 

In the typical analysis of observational data, it is necessary to discard samples that have been corrupted 
by experimental error, instrumental distortion, and/or radio frequency interference. This section demon- 
strates some of the PSRCHIVE algorithms that are available to assist in the automatic detection and 
excision of corrupted data. 

6.1 Excision of frequency channels that are known to be corrupted 

Although the observations of PSR J0437— 4715 used in this paper are not adversely affected by radio 
frequency interference, the data near the edges of the band are known to be corrupted by quantization 
distortions. 

Exercise: Use the psrsh command language interpreter and the zap edge command to assign zero 
weight to 15% of the total bandwidth (~ 19 frequency channels) on each edge of the band. Use the 
output option to write output data files with a new extension; e.g. 

cd $PSRCHIVE_DATA/mem 

psrsh - -e zz pulsar/*. ar cal/*.ar << EOD 

zap edge 0.15 

EOD 



In the above example, the single hyphen (-) command-line option instructs ps r sh to read the command 
script from the standard input. 

6.2 Automatic detection and excision of narrow-band interference 

Many types of radio frequency interference (RFI) are narrow band, such that the quality of an observa- 
tion can be significantly improved by discarding only a small number of corrupted frequency channels. 
The RFI environment at the telescope may be dynamic, such that it is not possible to select a fixed 
set of frequency channels to be excised at all times. To address this problem, PSRCHIVE implements 
an automatic frequency channel excision algorithm that is based on tolerance to differences between 
the observed spectrum and a version of the spectrum that has been smoothed by a running median. 
By default, the running median is computed using a window that is 21 frequency channels wide and 
all channels with total flux that differs from the median-smoothed spectrum by more than 4 times the 
standard deviation will be given zero weight. The standard deviation is defined recursively. That is, the 
algorithm works as follows 

1 . compute median-smoothed spectrum 

2. compute standard deviation, ignoring any zapped channels 

3. zap channels that differ from local median by more than tolerance 

4. if any channels were zapped, goto 2 

5. stop 

The above algorithm and its default parameters may not necessarily work in every situation, and it may 
require some experimentation to determine the parameters that detect the majority of RFI for a given 
telescope and instrument. 

Exercise: Use psrplot and the zap median pre-processing command to view an example of data 
corrupted by RFI that is not detected automatically by the default configuration of the algorithm de- 
scribed above. 

Cd $PSRCHIVE_DATA/zap/BPSR 

psrplot -p freq -j "zap median" example. ar 
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Fig. 1 Narrow-band radio frequency interference that is not automatically detected by the 
default configuration of the zap median algorithm. 

The plot produced by the above command is shown in Figure Q] Note that the signal above 1520 MHz 
(below channel index 10) has been filtered prior to digitization. Narrow -band RFI is evident just below 
1430 MHz and 1500 MHz. To understand why this seemingly-obvious RFI is not detected automatically, 
it is important to note that psrplot -p f req displays the pulsed flux as a function of pulse phase, 
whereas by default zap median works with the total flux summed over all pulse phases, which is 
plotted using 



psrplot -p psd example. ar 



and is shown in the top panel of Figure [2] Here, the channels that have been corrupted by RFI are not 
as obvious. To utilize a statistic that better characterizes pulsed flux, note that both the psd plot and the 
zap median algorithm can be configured to use any expression that is understood by psrstat. For 
example, 



psrplot -p psd -c ' exp= { $all : max-$all : min } ' example. ar 



produces the pulsed flux spectrum shown in the bottom panel of Figure [2] The same expression can be 
used to configure the zap median algorithm, as in the following psrsh script 

zap.psh 

#! /usr/bin/env psrsh 

# set the expression evaluated in each frequency channel 
zap median exp= { $all : max-$all : min } 

# execute the zap median algorithm 
zap median 

# zap frequency channels to 8 
zap chan 0-8 



This script can be passed to psrplot and used to process the data before plotting. 
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Fig. 2 Top panel: Spectrum formed by summing the values of all pulse phase bins for each 
frequency channel. Bottom panel: Pulsed spectrum formed by the difference between maxi- 
mum and minimum values of all pulse phase bins for each frequency channel. 



Cd $PSRCHIVE_DATA/zap/BPSR 

psrplot -p freq -J zap.psh example. ar 



Alternatively, the script can be made executable and run like a PSRCHIVE program with the standard 
output option to write the result to a file with a new extension. 

chmod a+x zap.psh 
./zap.psh -e zz example. ar 



6.3 Automatic detection and removal of impulsive interference 

Radio frequency interference (RFI) may also occur as broadband bursts of impulsive emission, such 
as lightning. When impulsive interference is persistent, it may not be possible to discard a subset of 
corrupted frequency channels or sub-integrations. To address this problem, PSRCHIVE implements an 
automatic impulsive interference mitigation algorithm that is based on tolerance to differences between 
the observed pulse profile and a version of the profile that has been smoothed by a running median. To 
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detect impulsive RFI of terrestrial origin, the pulse profile is first integrated over all frequency channels 
without correcting for interstellar dispersion. The running median is then computed using a window 
with a duty cycle of 2% and any phase bin with total flux that differs from the median-smoothed profile 
by more than 4 times the standard deviation is flagged for replacement. By default, the standard devi- 
ation is defined recursively in a manner similar to the algorithm used for excising corrupted frequency 
channels. The default recursive standard deviation estimator will fail in extreme cases of impulsive RFI, 
and in general it is better to enable the use of robust statistics (e.g. the median absolute deviation) as 
demonstrated in the exercise below. 

After flagging corrupted phase bins, the profiles in each frequency channel and polarization are 
corrected independently. A median-smoothed profile is computed and the values of flagged phase bins 
are set equal to the local median plus uncorrupted noise, defined as the difference between the observed 
profile and the local median of a randomly selected phase bin that has not been flagged as corrupted. 

Exercise: Use psrplot and the zap mow pre-processing command to view an example of data cor- 
rupted by impulsive RFI that is not detected automatically by the default configuration of the algorithm. 

Cd $PSRCHIVE_DATA/zap/DFB3 

psrplot -p freq+ -jp -j "zap mow" -1 subint=0- calibrator . ar 



The above command will produce four separate plots, one for each sub-integration, the first of which is 
shown in Figure [3] To enable the use of robust statistics and increase the median smoothing duty cycle 
from 2% to 10%, create the following psr sh script 

mow.psh 

#! /usr/bin/env psrsh 

# use robust statistics (median absolute deviation) 
zap mow robust 

# set the median-smoothing window to 10% of the pulse profile 
zap mow window=0 . 1 

# execute the zap mow algorithm 
zap mow 



and pass this script to psrplot to process the data before plotting. 

Cd $PSRCHIVE_DATA/zap/DFB3 

psrplot -p freq+ -jp -J mow.psh -1 subint=0- calibrator . ar 



The first plot produced by running the above command is shown in Figure |4] 
6.4 Interactive excision with p s r z ap 

In some cases, automated methods of detecting corrupted data are insufficient and it is necessary to 
perform the task manually. Large data files with multiple sub -integrations and frequency channels are 
best visualized using a dynamic spectrum in which the colour of each pixel in a two-dimensional im- 
age spanned by time and frequency is determined by a value computed from the pulse profile at that 
coordinate. The psrzap program is an interactive tool for excising corrupted data using the dynamic 
spectrum. Run psrzap -h for a list of the keyboard and mouse interactive commands, then 

Cd $PSRCHIVE_DATA/zap/GUPPI 

psrzap guppi_552 4 5_190 9-37 4 4_0033_0001 . rf 
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Fig. 3 Impulsive radio frequency interference that is not automatically detected by the default 
configuration of the zap mow algorithm. 
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Fig. 4 Impulsive radio frequency interference is detected by the zap mow algorithm after 
enabling the use of robust statistics and increasing the median smoothing duty cycle to 10%. 



After loading the data file, two plot windows will open. The main window plots the dynamic noise 
spectrum, which by default is defined as the variance in each pulse profile as a function of sub-integration 
(x-axis) and frequency (y-axis). The secondary plot window is divided into three panels: the top panel 
displays the pulse profile after integration over all time sub-integrations and frequency channels, the 
middle panel displays the phase-versus-frequency image after integration over all time sub-integrations, 
and the bottom panel displays the phase-versus-time image after integration over all frequency channels. 
These diagnostic plots are updated by pressing d on the keyboard. 

As with psrplot -p psd and zap median, the quantity that is plotted as a function of time 
and frequency may be specified using a mathematical expression that is understood by psrstat; e.g. 



psrzap -E ' { $all :max-$all :min} ' guppi_55245_1909-3744_0033_0001 . rf 
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There are three modes for selecting ranges of data to view or excise: 

1 . time (t on keyboard) selects an entire sub-integration (column) 

2. frequency (f on keyboard) selects an entire frequency channel (row) 

3. both (b on keyboard) selects a rectangular region 

The line(s) passing through the cursor indicate the current selection mode. To zoom in on a desired 
range, click the left mouse button at the start and end positions of the range. To excise a desired range, 
click the left mouse button at the start, and the right mouse button at the end of the range. Simply right 
clicking will excise the single column, row, or pixel under the mouse (depending on the selection mode). 

RFI typically appears as bright spots in the dynamic noise spectrum. Excise corrupted data until you 
are satisfied and save the result by pressing s on the keyboard. Press w to generate a psrsh command 
script that reproduces the results of the interactive excision session; the script is saved as a text file with 
a filename created by appending the . psh extension to the filename of the input archive. This script can 
be integrated into an automated pipeline that reprocesses the original data from scratch. Use psr stat 
to compare the S/N of the data before and after RFI excision. 



psrstat - jTFp -c snr guppi_55245_1909-3744_0033_0001 . rf * 



7 POLARIMETRIC CALIBRATION 

Polarization measurements provide additional insight into the physics of both the emission and propa- 
gation of electromagnetic radiation; e.g. measurements of Fara day rotation in th e interstellar medium 
yield constraints on the structure of the Galactic magnetic field (lHan et all 120061) and estimates of the 
position angle of the linearly polariz ed flux indicate that the pulsar spin axis may be aligned with its 
space velocity (Ijohnston et all 120051) . It has also been demonstrated that accurate polarimetry and ar- 
rival time estimation using all four Stokes param eters can signficantly improve both the accuracy and 
precision of pulsar timing data (Ivan Stratenll2006l) . 

The processes of reception and detection introduce instrumental artifacts that must be calibrated 
before meaningful interpretations of experimental data can be made. A first-order approximation to cal- 
ibration can be performed using observations of a noise diode that is coupled to the receptors, from 
which the complex gains of the instrumental response as a function of frequency are derived. This ap- 
proximation to calibration is based on the ideal feed assumption that the Jones matrix in each frequency 
channel has the form 



where zq and z\ are the complex gains. The absolute phase of the Jones matrix is lost during detection, 
and the matrix may be parameterized using a polar decomposition described by the absolute gain G, 
differential gain 7, and differential phase </>. 

7.1 Display calibrator parameters with p s rp 1 ot 

In PSRCHIVE, calibrator observations of the noise diode have type=PolnCal (as returned by 
psredit); the noise diode is typically driven by a square wave with a 50% duty cycle, as shown in 
Figures |3] and H 

Exercise: Plot the polar decomposition of the ideal feed as a function of frequency. 



cd $PSRCHIVE_DATA/mem/cal 
psrplot -p calm *.zz 




(1) 
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Fig. 5 Observations of the Parkes 21 -cm Multibeam receiver noise source. The total intensity 
from a single 500kHz channel was integrated for approximately 80s. In the top panel, the 
telescope was pointed at 3C 218 (Hydra A); the mean on-pulse power (green) is used to esti- 
mate H on and the mean off-pulse power (red) is used to estimate L on . In the bottom panel, 
the telescope was pointed 2 deg north; the mean on-pulse power is used to estimate H Q s and 
the mean off -pulse power is used to estimate L s . 



7.2 Prepare flux calibrator data with f luxcal 

Although not immediately necessary for pulsar timing, it is also useful to perform absolute flux 
calibration; for example, well-c alibrated estimates o f flux density may be used in long-term studies 
of refractive scintillation (e.g. iRickett et all Il984l) . Absolute flux calibration is performed using 
observations of a standard candle, an astronomical source with a well-determined reference flux 
density and spectral index that applies over a broad range of radio frequencies. Two sets of observations 
are made: the noise diode is driven while the telescope is pointed at 1) the standard candle and 2) 
a nearby patch of sky that is assumed to be empty. Example plots derived from the sample data are 
shown in Figure [5] Note that the integration lengths for the on- and off-source observations need not 
necessarily be equal (as long as the data represent mean flux densities). The absolute gains also need 
not be equal. Given 
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Hon — gon(S sys + T + Co) H g = g s(S sys + Co) 

ion = . < 7on(S'sys + Tq) L a ff = g ffS S ys 

where g on and <7 Q ff are the unknown absolute gains of the instrument while pointing on and off the 
standard candle, S sys is the unknown system equivalent flux density, So is the known flux density of the 
standard candle, and Co is the unknown flux density of the receiver noise source. Then, 

/on = ~J 1 = r, (2) 

^on >3sys T tJQ 



and 



/off - j— - 1 - -z— (3) 

-<^off ^sys 



and 

±-± = ?l (4) 

/on /off C 



Equation [4] is solved for Co, then Equation [3] is solved for 5 sys . In PSRCHIVE, flux calibration 
observations have a type attribute equal to FluxCalOn or FluxCalOf f . 

Exercise: Start by creating a calibrator database. 

cd $PSRCHIVE_DATA/mem/cal 
pac -w -u zz 



This creates a file called database . txt. Then run 



fluxcal -f -d database.txt 



This will produce a file named n2003201035947 . fluxcal and update database .txt with a 
new entry for this file. Use psrplot -p calm to plot the derived estimates of S sys and Sq as a 
function of radio frequency). 

7.3 Correct the receiver parameters 

A large number of assumptions are built into the design of an instrument, ranging from the sign of the 
complex argument in exp(±iwi) to the handedness of circular polar ization. Over time, a numb er of 
inconsitent conventions have been utilised by various authors; e.g. see lEverett & Weisbergl (1200 ll) for a 
thorough review of contrad icting definitions o f the p osition angle of the linearly polarized flux. 

To address this issue, Ivan Straten et al.l (1201 Ol) define the PSR/IEEE convention that is used by 
PSRCHIVE and include a table of parameters that can be used to describe the differences between an 
instrumental design and the PSR/IEEE convention. For the Parkes 21 -cm Multibeam receiver, it is nec- 
essary to set the symmetry angle to — n/2, which can be done with the following command 



cd $PSRCHIVE_DATA/mem 

psredit -c rcvr:sa=-90 -m pulsar/*. ar cal/*.ar 



7.4 Calibrate using the ideal feed assumption 

To perform the first-order approximation to calibration based on the ideal feed assumption, run 
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cd $PSRCHIVE_DATA/mem/pulsar 
pac -d . . /cal/database . txt *.zz 



For each input file, a new output file will be written with a filename created by appending the extension 
. calib to the input filename. If the receiver were ideal, the first-order approximation to calibration 
would have eliminated the variation of the Stokes parameters as a function of parallactic angle. Use 
psrplot to test this expectation. 



psrplot -ps -jF *. calib 



The Stokes parameters still vary as a function of time because the ideal feed assumption does not apply 
to the Parkes 21-cm Multibeam receiver. Create a sub-directory, e.g. ideal/ and move the newly 
calibrated data to this sub-directory (otherwise, they will be over-written in the next step). 



mkdir ideal 

mv *. calib ideal/ 



7.5 Measure the cross-coupling parameters using pern 

To accurately calibrate these data, the cross-coupling terms (off-diagonal components of the Jones ma- 
trix) must be estimated, which can be done by modeling the variation of the Stokes parameters as a func- 
tion of time. The process of performing the least-squar es fit is called Mea surement Equation Modeling 
(MEM); the PSRCHIVE implementation is described in Ivan Stratenl (120041) and more fully documented 
onlinq3- Use psradd to combine the archives calibrated using the ideal feed assumption into a single 
archive 

psradd -T -o calib. TT ideal/* . calib 



The resulting archive will be used as input to pem, from which it will choose the best phase bins to use 
as constraints and derive the first guess for the polarization of the source. 

pem -d .. /cal/database . txt -s -c calib. TT *.zz 



While running, pem outputs messages about the quality of the least-squares fits, which are performed 
independently in each frequency channel. On a multi-processor machine, multiple channels may be 
solved simultaneously by using the -t <nthread> command-line option, where <nthread> is 
the number of processing threads to run in parallel. When pem finishes, it produces an output file 
pem . f it s that contains the MEM solution; the model parameters may be plotted using 



psrplot -p calm pem. fits 



Compared to the solution derived using the ideal feed assumption, three new parameters have been 
added to the model of the receiver: 9\ describes the non-orthogonality of the feed receptors (the linearly 
polarized receptors should be oriented at and 90 degrees) and are the ellipticities of the receptors, 
which should be in an ideal feed with linearly polarized receptors. The mean value of ~ 5 degrees 
corresponds to roughly 15% mixing between linear and circular polarizations (Stokes Q and V). It is 
not possible to determine the absolute rotation of the receptors about the line of sight, 9q, without an 
external reference; therefore, only the non-orthogonality is measured. 



http://psrchive.sourceforge.net/manuals/pcm 
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The solution output by pcm also includes estimates of the Stokes parameters of the noise diode, 
which is no longer assumed to illuminate both receptors equally and in phase. This information enables 
the calibrator so lution derived from one data set to be applied to observations of another source, as in 
IOrdetal.l(l2004l) . 

7.6 Calibrate using the MEM solution 

Move the output file pcm. fits to the cal/ sub-directory, change to this directory, and recreate the 
calibrator database 

mv pcm. fits ../cal/ 
cd . . /cal 

pac -w -u zz -u fits -u fluxcal 



Confirm that pcm . fit s has been added to database . txt, then run 
cd . . /pulsar 

pac -d . . /cal/database . txt -S *.zz 



Plot the Stokes polarization profile (integrated over the entire band) in each calibrated data file output 
by pac to confirm that the Stokes parameters no longer vary with time. 

8 ARRIVAL TIME ESTIMATION 

In this section, arrival time estimates are derived using a previously created standard (template) profile. 
This high S/N standard profile was formed by integrating data from ~ 42 hours of observations. 

8.1 Prepare the standard profile 

The standard profile is located in $PSRCHIVE_DATA/mem/ std/ standard . ar. Plot the phase-vs- 
frequency image of the total intensity and compare this image with that of nonspc . ar in the same 
directory. The file nonspc . ar was formed from data that were not corrected for scattered power. 
Although standard . ar was corrected, there are still residual artefacts in the edges of the band. Use 
p s r s h to give zero weight to the affected frequency channels then integrate over all frequency channels. 
For best results, excise the same frequency channels that were excised from the data in Section |6*T1 e.g. 

cd $PSRCHIVE_DATA/mem/std 

psrsh - -e FF standard. ar << EOD 

zap edge 0.15 

f scrunch 

EOD 



If the GNU Scientific LibrarvPI is installed and detected during the configuration of the PSRCHIVE 
software, then it is possible to use the wavelet smoothing algorithm implemented by psrsmooth to 
create a "noise-free" template profile. 



psrsmooth -W -t UD8 standard. FF 



This will produce a file called standard . FF . sm. By default, psrsmooth a pplies the translation- 
invariant wavelet denoising algorithm described by ICoifman & Donohol d 1995b . The profile data are 



http://www.gnu.org/software/gsl 
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first transformed into the wavelet domain, where a noise level is estimated from the data. Based on the 
measured noise level, a threshold is calculated, and all wavelet coefficients with absolute value below 
the threshold level are set to zero. The data are then transformed back into the profile domain, resulting 
in a smoothed profile. Use the crop attribute of the flux plot to zoom in on the low amplitude flux 
near the off -pulse baseline and compare the standard profile with its smoothed version; e.g. 

psrplot -pD -jp -c crop=0.01 -N 1x2 standard. FF standard . FF . sm 



8.2 Estimate arrival times using pat 

In this section, two different methods of arr ival time estimation are compared: scalar template matching 
using only the total inten sity dTaylorl["l992l) and matrix template matching using all four Stokes param- 
eters (Ivan StratenL 120061) . Comparison is also made between the results derived using the two different 
template profiles: the smoothed and not smoothed versions of standard . ar. Finally, there are three 
different data sets: the uncalibrated data, the data calibrated using the ideal feed assumption, and the 
data calibrated using the MEM solution derived with pcm. In total, there are 12 different combinations 
of arrival time estimation algorithm, template profile, and observational data. Experiment with these 
combinations to find the arrival times with the lowest residual standard deviation. To experiment, run 
pat in either scalar template matching mode; e.g. 

cd $PSRCHIVE_DATA/mem/pulsar 

pat -F -s . . /std/standard . FF *.zz > uncal_unsmooth_stm . tim 



or matrix template matching mode; e.g. 

pat -Fpc -s .. /std/standard . FF . sm *.calib > cal_smooth_mtm . tim 



Run tempo2 to evaluate the arrival times; e.g. 

tempo2 -f .. /pulsar . par uncal_unsmooth_stm . tim 



Search the output of tempo 2 for lines like 



RMS pre-fit 


residual = 


0.11 (us), RMS 


post- 


-fit residual = 0.11 (us) 


Fit Chisq = 


359.7 


Chisqr/ nf ree 


= 359 


74/95 = 3.78671 



and make note of both the RMS post-fit residual and Chisqr/nfree in each case tested. 



8.3 Correct arrival time estimation bias using psrpca 

As the mean flux density of a source of noise approaches the system equivalent flux den sity, the statis- 
tics of the noise intrinsic to the source (or self-noise) can no longer be neglected (e.g. iGwinnl 12001 1 
van StratenL l2009t lG winn & Johnsonll201 lb . This i s particularly t rue when the signal is heavily modu- 
lated, as is typically the case with pulsar emission (Ric kets | 19751). which can be described as stochas- 
tic wide-band impulse modulated self-noise (SWIMS: lOslowski et all 1201 ll) . Pulsed self-noise is het- 
eroscedastic and, when the timescale of impulsive modulation is longer than the sampling interval re- 
quired to resolve the mean pulse profile, SWIMS is correlated. Correlated and heteroscedastic noise 
violates the basic premises of least-squares estimation and introduces pulse arrival time measurement 
bias. This bias may be corrected using psrpc aFl which performs a principle component analysis of 
the observed pulse profile shape fluctuations and a multiple regression analysis in which the post-fit 



The psrpca program will be compiled only if the GNU Scientific Library is installed 
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arrival time residual is the dependent variable and the most sig nificant principal com ponents are the 
indepe ndent variables. The methodology is described in detail by iDemores 3 d2007l) and Oslowski et al. 
(1201 lh . 

SWIMS is best characterised using a large quantity of data; at the very least, the number of ob- 
servations must exceed the number of phase bins used to resolve the mean pulse profile. This con- 
straint is satisfied by the example data: ~ 9000 pulse profiles spanning ~ 10 days of observations of 
PSR J0437 —47 1 5 made at 20 cm with the Parkes 64 m radio telescope between 1 9 and 27 July 2003 . To 
measure and remove the arrival time bias introduced by SWIMS, first produce the arrival time estimates 



Cd $PSRCHIVE_DATA/pca 
Is -1 *.ar > files. Is 

pat -s . . /mem/std/standard . FF -f tempo2 - 


-M files 


.Is > psrpca.tim 


then compute the post-fit arrival time residuals 


tempo2 -output general2 -s '{sat} {post} {err} 
-f .. /mem/pulsar . par psrpca.tim | grep SWIMS \ 
awk '{print $1,$2,$3}' > resid.dat 


SWIMS\n' \ 


and perform the principal component and multiple regression 


analyses 




psrpca -s . . /mem/std/stanard . FF -r resid 


dat -M 


files . Is 



When running psrpca, it is important to ensure that the arrival time residuals and input data files 
are provided in exactly the same order. The above commands ensure this by first creating a file listing 
named files . Is, which is passed to both pat and psrpca. A variety of diagnostic output files are 
produced by psrpca, each with a name that starts with the prefix psrpca; this prefix can be chosen 
by using the -p prefix option. The diagnostic output files are: 

- psrpca_dif f s . ar - an archive containing the differences between the standard template and 
observed pulse profiles; these data are used to construct the covariance matrix; 

- psrpca_covariance . dat - a plain text file containing the covariance matrix; 

- psrpca_evals . dat - a plain text file containing the eigenvalues derived from the covariance 
matrix in a format that is easily inspected using GNUPLOT; 

- psrpca_evecs . ar - a PSRFITS archive containing the eigenvectors; 

- psrpca_decomposition.dat - a plain text file containing the decompositions of the ob- 
served pulse profiles onto the measured eigenvectors; 

- psrpca_beta_zero.dat and psrpca_beta_vector_used.dat - plain text files con- 
taining the regression coefficients used to remove the bias in arrival time estimates; and 

- psrpca_residuals . dat - a plain text file containing the bias-corrected arrival time residuals. 

These output files enable easy inspection of the principal component analysis results. The plain-text 
output file named psrpca_residuals . dat contains four columns: 

1. MJD - the site arrival time; 

2. biased_residual - the arrival time residual produced by tempo2 (expressed in fj,s) 

3. corrected_residual - the bias-corrected arrival time residual 

4. error - the estimated measurement error produced by pat 



This file can be used to compare the biased and corrected arrival time residuals. The provided script 
rms . sh may be used to compare the standard deviation of the two sets of arrival time residuals; e.g. 



. / rms 


sh 


psrpca_residuals 


dat 


2 4 


. / rms 


sh 


psrpca_residuals 


dat 


3 4 
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The first command produces the weighted standard deviation and a measure of goodness of fit for the 
biased arrival time residuals; the latter computes the same for the bias-corrected residuals. The biased 
and corrected residuals may also be inspected using GNUPLOT. 

The primary output of the psrpca program is a file named psrpca_std . ar. This PSRFITS 
archive contains a copy of the standard template profile that was provided to psrpca with an addi- 
tional extension that contains the principal component eigenvectors and the multiple regression coeffi- 
cients, pat can use the information in this extension to correct the bias in output arrival time estimates. 
With this functionality, the bias predictor derived from one epoch may be applied to observations made 
at other epochs and the bias-corrected arrival time estimates may be provided as input to TEMPO or 
TEMP02, thereby yielding improved physical parameter estimates. 

9 DISCUSSION 

This paper presents a scientific workflow for high-precision timing experiments and demonstrates some 
of the general-purpose tools applicable to a wider variety of pulsar studies. In addition to the programs 
described here, PSRCHIVE development continues on a number of novel anal ysis tools. 

F or example, the psrmode l prograrrPI fits the rotating vector model (iRadhakrishnan & Cookei 
1 19691 lEverett & Weisbergl 12001 1) to observed polarisation data using a statistically robust algorithm. 
Rather than perform a one-dimensional fit to the real-valued position angle as a function of pulse phase, 
psrmodel performs a two-dimensional fit directly to the Stokes Q and Stokes U profiles by treating 
them as the real and imaginary components of a complex number. For each complex number, the phase is 
given by the position angle predicted by the rotating vector model and the magnitude (linearly polarized 
flux) is modeled as a free parameter. This approach has a number of advantages over directly modeling 
the position angle: 1) because Stokes Q and U errors are normally distributed, low S/N data may be 
included in the fit without biasing the result; 2) Stokes Q and U are not cyclic; and 3) orthogonal mode 
transitions are trivially modeled by negative values of the complex magnitudes. 

Also currently under development, the psrspa programrjcan search for significant single pulses 
and derive a wide variety of statistic al quantities from sing le-pulse data, such as the phase-resolved 
histogram of the position angle (e. g. IStinebring et aU 1 1984)) and the two-dimens ional distribution of 
the polarization vector orientation (lEdwards & Stappersi 120041 iMcKinnorl 1 20091) . This program was 
recently used to a nalyse the phase d istribution and width of single pulses from the radio magnetar, 
PSR J1622-4950 (iLevin et all 12012)) . 

A number of the algorithms implemented by PSRCHIVE consist of a single, independent profile 
transformation that is performed sequentially by looping over all sub-integrations and frequency chan- 
nels. On a multiprocessor architecture, such transformations could be readily executed in parallel; fur- 
thermore, identification of such parallelism also provides the oppor tunity to apply loop transformations 
that i mprove data locality and conserve memory bandwidth (e.g. lAbu-Sufahl 1 19791: iMcKinlev et all 
119961) . 

Though written in C++, it is possible to access a large fraction of PSRCHIVE functionality via 
the Python programming languagq£j This provides an alternative to the PSRCHIVE applications and 
psrsh command language interpeter for both high-level scripting of PSRCHIVE functionality and 
interactive or non-standard data manipulation. The Python interface works by providing direct access 
to the core C++ class library on which PSRCHIVE is built. Additional information about installing 
and using the Python interface is available in the online documentation^ Developers with an interest 
in extending, refining, or optimising the PSRCHIVE software are encouraged to contact the project 
administrators and refer to the extensive online documentation hosted by SOURCEFORGrJ^) 



http://psrchive.sourceforge.net/manuals/psnTiodel 
http://psrchive.sourceforge.net/manuals/psrspa 
http ://w ww. python . org 

http://psrchive.sourceforge.net/manuals/python 
http://psrchive.sourceforge.net/devel 
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Appendix A: OBTAINING DATA FROM THE CSIRO DATA ACCESS PORTAL 

To obtain a copy of the observational data used in this paper, visit the CSIRO Data Access Portal and, 
using the Pulsar Search tool, enter the following information 

- Source Name: J0437-47 15 

- Project ID: PI 40 

- Observation Date (dd/mm/yyyy): 19/07/2003 to 27/07/2003 

Using the check boxes in the column on the left, refine the search results to include only "raw" obser- 
vations (i.e. not preprocessed) with a frequency of 1341 MHz. This should yield 67 results, the first of 
which has the filename n2003-07-19-18:08:01.rf,anda total download size of 1 5 . 3 GB . 

The data stored on the CSIRO Data Access Portal have more time and frequency resolution than is 
required for the purposes of the demonstrations presented in this paper. Processed versions of these data 
with lower resolution are also available for download from Swinburne University of Technology at 

\protect \vrule width Opt \ protect \ href { http : // astronomy . swm .edu.au/pulsar/ dat 

Here, there are three files 

- psrchive_mem . tgz (112.6 MB) contains 5-minute integrations with full frequency resolution 
from the first day of observations (19 July 2003) - these data are for use in the exercises presented 
in Sections[2]through|8j 

- psrchive_zap . tgz (299.4 MB) contains three files that are not available via the CSIRO Data 
Access Portal - these data are for use in the exercises presented in Section|6l and 

- psrchive_pca . tgz (132.4 MB) contains ~ 17-second integrations with no frequency resolu- 
tion spanning the nine days of observations (from 19 to 27 July 2003) - these data are for use in the 
exercises presented in Section liOl 

These files can be unpacked with commands such as 



gunzip -c psrchive_mem . tgz | tar xf - 



The examples presented throughout this paper assume that the above three files have been downloaded 
from Swinburne University of Technology and the full path to the directory into which they have been 
unpacked is recorded using the $PSRCHIVE_DATA shell environment variable. 
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